Binary black-hole evolutions of excision and puncture data 
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We present a new numerical code developed for the evolution of binary black-hole spacetimes 
using different initial data and evolution techniques. The code is demonstrated to produce state-of- 
the-art simulations of orbiting and inspiralling black-hole binaries with convergent waveforms. We 
also present the first detailed study of the dependence of gravitational waveforms resulting from 
three-dimensional evolutions of different types of initial data. For this purpose we compare the 
waveforms generated by head-on collisions of superposed Kerr-Schild, Misner and Brill-Lindquist 
data over a wide range of initial separations. 
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I. INTRODUCTION 

In the course of the last two years, the research area 
of gravitational wave physics has entered a very excit- 
ing era. On the experimental side, the first generation 
of ground-based Gravitational Wave detectors, LIGO, 
GEO600, TAMA300 and VIRGO, are performing obser- 
vation runs at, and even beyond design sensitivity [1-5]. 
At the same time, the simulation of the most promising 
sources of gravitational waves, the inspiral of compact 
binary systems, has made enormous progress. While ap- 
proximate studies based on the Post-Newtonian approach 
have been able for some time to accurately simulate the 
earlier stages of inspiralling binary systems [6-11], recent 
developments in numerical relativity have made possible 
the simulation of the highly relativistic final stages of the 
inspiral and merger of compact binaries in the framework 
of fully nonlinear general relativity. 

For a long time such simulations have been troubled 
by stability problems which caused evolutions to termi- 
nate after times relatively short compared with the dy- 
namic timescales of the problems under investigation. It 
is becoming increasingly clear now, however, that these 
problems have been successfully overcome by a combina- 
tion of modified formulations of the Einstein equations 
[12-15], suitable gauge conditions (see e.g. [16-19]) and 
improved techniques for the treatment of the singularities 
inherent to black-hole spacetimes. 

Using such modern techniques, Briigmann et. al. [17] 
obtained the first simulation of a complete orbit of a 
black-hole binary in the framework of the Baumgarte- 
Shapiro-Shibata-Nakamura (BSSN) formulation [12, 13], 
using puncture data [20] and corotating coordinates. 
More recently, their results have been confirmed by an 
improved study [21]. The first waveforms generated in 
the inspiral and coalescence of black holes have been 
presented by Pretorius [22, 23] who uses a generalized 
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harmonic formulation of the Einstein equations com- 
bined with special numerical techniques such as black- 
hole excision, spatial compactification and implicit fi- 
nite differencing. The latest development, simultane- 
ously discovered by Campanelli et. al. and Baker et. al. 
[18, 19, 24, 25], is based on the evolution of black-hole 
data of puncture type using special gauge conditions ac- 
commodating the motion of the punctures across the 
computational domain. For this reason these simulations 
are commonly referred to as moving punctures. More re- 
cently, this technique has facilitated the investigation of 
various aspects of the binary-black-hole coalescence, such 
as the radiation of linear momentum by systems of un- 
equal masses and/or spins [26-34], the impact on the 
waveforms and merger dynamics of nonvanishing spins 
[35, 36] and analysis of the waveforms in the framework 
of post-Newtonian inspiral and black- hole ring-down [37- 
42]. 

As in the case of black holes, simulations lasting for 
several orbits have also been obtained for neutron star bi- 
naries by several groups [43-46] . In more recent develop- 
ments the focus is switching to the refinement of the mat- 
ter models, as for example by the inclusion of magneto- 
hydrodynamic effects (see e.g. [47]). In our work, how- 
ever, we focus on black-hole systems, and will therefore 
exclusively study vacuum spacetimes. 

In spite of the dramatic progress in numerical simula- 
tions of black-hole binaries, there remain important ques- 
tions to be answered, in particular with regard to the use 
of the resulting waveforms in the ongoing effort to detect 
and physically interprete gravitational- wave signals. In 
particular it will be important to establish the accuracy 
of the numerically calculated waveforms and the consis- 
tency of these results with regard to the use of different 
types of binary-black-hole initial data and the evolution 
techniques used in the codes. First steps in this direction 
have been undertaken with regard to the use of evolution 
techniques and separation parameters of a given initial- 
data type. In Ref. [48], the impact of black-hole exci- 
sion was studied in the case of head-on collisions of Brill- 
Lindquist data. The results with and without excision 
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yielded good agreement in that study. A comparison be- 
tween plunge waveforms obtained from moving-puncture 
evolutions with those resulting from Lazarus calculations 
[49] has been presented in Ref. [19]. Furthermore, the 
waveforms resulting from inspiralling black holes of punc- 
ture type starting from different separations have been 
found to show excellent agreement in Refs. [27, 35]. A 
comparison of waveforms obtained from evolving concep- 
tually different types of initial data has been presented in 
[50]. This study is inhibited, however, by the difficulties 
in starting the simulations from comparable initial con- 
figurations as is demonstrated by the nonvanishing spin 
in one of the two data sets considered in this work. 

The purpose of this paper is two-fold. First, we present 
a new numerical code which has been designed to ac- 
commodate different types of initial data, formulations of 
the Einstein equations as well as singularity treatment. 
We demonstrate that the code is capable of producing 
state-of-the-art simulations of inspiralling black-hole bi- 
naries and extract convergent waveforms. Second, we 
use the code to further progress in the comparison of 
different initial configurations by comparing black-hole 
head-on collisions obtained from different types of initial 
data and using different evolution techniques. Specifi- 
cally, we compare the results obtained from superposed 
Kerr-Schild data evolved in the framework of black-hole 
excision and algebraic gauge conditions with those ob- 
tained from evolving Brill-Lindquist as well as Misner 
data in the framework of the moving-puncture method. 

This paper is structured as follows. We begin with 
a detailed presentation of the code in Sec. II. Next, we 
benchmark the code in Sec. Ill by simulating the inspiral 
and merger of an orbiting black-hole binary comparable 
to those studied in the recent literature. The compar- 
ison of head-on collisions obtained with Brill-Lindquist, 
Misner and Kerr-Schild data is given in Sec. IV and we 
conclude with a discussion of our findings in Sec. V. De- 
tails on the exact version of the BSSN equations used 
for this work, the analytic solution of a boosted black 
hole in Kerr-Schild coordinates, extraction of gravita- 
tional waves and the performance of the code are pre- 
sented in Appendices A-D. Throughout this work we set 
G = c = 1 and use greek indices for spacetime compo- 
nents 0...3 and latin indices for spatial components 1...3. 



II. COMPUTATIONAL FRAMEWORK 

The simulations presented in this work have been ob- 
tained with the newly developed Lean code. This code 
has been inspired partly by the Maya^ code [51-53], and 
partly by the most recent developments in the simula- 



Throughout this work with the Maya code we refer to the version 
used in Rcf. [51] which is not to be confused with the new code 
of the same name described in Ref. [26] . 



tion of black- hole data of puncture type [24, 25]. It is 
based on the Cactus computational toolkit [54], used for 
parallelization and data input/output. Mesh refinement 
is provided by Carpet [55, 56], puncture initial data 
by the TwoPunctures thorn [57] and horizon finding 
by AHFinderDirect [58, 59]. The code achieves dy- 
namic mesh refinement by stcxTing in accordance with 
the black-hole motion the regridding option inherent to 
the Carpet package. While the Lean code has been in- 
spired by Maya, it has been written entirely from scratch 
and various new features have been added. These are 
fourth-order discretization of the spatial derivatives, the 
evolution of the BSSN equations using the x version (see 
below), additional gauge conditions, time integration us- 
ing the fourth-order Runge-Kutta (RK) scheme, dynamic 
mesh refinement that allows for multiple refinement com- 
ponents to follow the black-hole motion and merge into 
single components and additional initial-data options in- 
cluding puncture data using the TwoPuncture thorn 
[57] and Misner data. Furthermore, the different orga- 
nization of the code has lead to about five times faster 
evolutions and a reduction by about a third in memory 
requirements compared with Maya for a given config- 
uration. Details on the code's performance for the or- 
bital simulations and head-on collisions are provided in 
App. D. The key feature of the code for the compari- 
son presented below is the incorporation in the frame- 
work of mesh refinement of both, dynamic black-hole 
excision and the moving-puncture technique used with 
enormous success in evolutions of conformally flat initial 
data. These features as well as other aspects of the Lean 
code are described in more detail in the remainder of this 
section. 



A. Formulation of the Einstein equations 

Most of the numerical work in three spatial dimensions 
has been performed inside the framework of the canon- 
ical "3+1" spacetime decomposition of Arnowitt, Descr 
and Misner (ADM) [60] (see also [61] for a detailed dis- 
cussion). In the notation of [61], the geometry is de- 
scribed in terms of the three-dimensional metric 7y and 
the extrinsic curvature JsT^, as well as four gauge func- 
tions a and /3* which represent the coordinate freedom 
of general relativity. The Einstein field equations result 
in six evolution equations each for and Kij as well as 
four constraint equations, namely the Hamiltonian and 
momentum constraints. These equations are commonly 
referred to as the ADM equations. 

While these equations have been at the heart of most 
numerical codes for a long time, the ensuing stability 
problems have lead to the use of various alternative for- 
mulations of the Einstein equations, most of them mod- 
ifications of the ADM equations. The most popular and 
successful of these modified schemes is now known as the 
BSSN system [12, 13] and has been implemented in the 
Lean code. While the code also allows evolutions us- 
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ing the Nagy-Ortiz-Reula (NOR) [62] or the generalized 
harmonic formulation [14, 15], we have not yet managed 
to achieve long-term stable simulations using these sys- 
tems. Therefore, all simulations presented in this work 
have been obtained with the BSSN system. 

The BSSN formulation results from applying the fol- 
lowing modifications to the original ADM equations: 
First, a split of the extrinsic curvature into a trace- free 
part Aij and the trace K, second, a conformal rescal- 
ing of the three-metric and the extrinsic curvature and, 
third, the introduction of contracted Christoffel symbols 
as separate variables r\ One thus arrives at a description 
of the spacetime in terms of the variables 

(/>=j^ln7, 7jj- = e~"*'^7y. 



•pi — ;r,mnpi 



(1) 



as well as the gauge functions a and Here, 7 denotes 
det 7jj and the definition of cj) implies that 7 = det 7^ = 
1. Alternatively to this choice of variables the Lean code 
also allows for evolutions using the variable 



X 



(2) 



as introduced in Ref . [24] . The complete evolution equa- 
tions for both sets of variables are listed in Appendix A. 

We refer to the two resulting systems given by Eqs. (Al)- 
(A5) and Eqs. (Al), (A4), (A6)-(A8) as the cj) and the 
X version of the BSSN equations respectively in the re- 
mainder of this work. 

In addition to evolving the BSSN variables according 
to either of these systems, we enforce after each update 
of the variables the condition A% = 0, which is a conse- 
quence of the definition of Aij in Eq. (1). We find this 
step to be crucial for the stability of our simulations. 
Other modifications to the BSSN equations have been 
suggested in the literature [see e.g. [63]]. We have ex- 
perimented with several of these, but not observed any 
further improvements of the performance of the code. In 
particular we do not find it necessary to enforce the con- 
dition det "fij = 1 or to replace the variable in terms 
of the Christoffel symbols at any stage of the evolution. 



B. Initial data 

One main purpose of this paper is to provide a detailed 
comparison of binary-black-hole collisions obtained with 
different initial-data types. We now describe the different 
initial data available inside the code. Specifically, these 
are puncture, Misner and superposed Kerr-Schild data. 

The starting point for binary-black-hole data of the 
puncture type is the Schwarzschild solution in isotropic 
coordinates, where the spacetime curvature is captured 
entirely within the conformal factor tp = e* = (14-^). 
In the case of time symmetry, these conformally fiat data 



have been shown to generalize to an arbitrary number of 
black holes by merely adding the individual quotients in 
the conformal factor [64, 65] 



^ = 1 + 



2|r - fi\ 



(3) 



where the index i labels the individual black holes. This 
time-symmetric initial configuration of multiple black 
holes is known as Brill-Lindquist data. As a further gen- 
eralization of these data, spin and momentum can be 
incorporated in the form of a nonvanishing extrinsic cur- 
vature [66]. Finally, Brandt and Briigmann [20] have 
transformed this type of data into a form substantially 
more convenient for the use in numerical simulations by 
applying a compactification to the internal asymptoti- 
cally fiat regions of the holes (see their paper for exis- 
tence and uniqueness of the solutions for the Hamilto- 
nian constraint). These data are commonly referred to 
as punctures and have been widely used in numerical 
simulations. 

Inside the Lean code, initial data of Brill-Lindquist 
type are implemented analytically using Eq. (3). More 
general classes of puncture data are made available via 
the TwoPuNCTURES thorn of Ansorg et. al. [57], which 
solves the Hamiltonian constraint using spectral methods 
combined with transformations to a coordinate system 
specially adapted to the structure of the binary-black- 
hole spacetime [see [57] for details]. 

The second class of initial data we study in this work 
are the axisymmetric Misner data [67] which represent 
a conformally fiat spacetime containing two nonspinning 
equal-mass black holes at the moment of time symmetry 
[Kij =0). In Cartesian coordinates the three-metric jij 
for this configuration can be written as 



where the conformal factor is given by 



(4) 
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v/.t2 + ?y2 + (2; _ z„y 
Zn = coth n/i. 



(5) 

(6) 



and is a free parameter determining the initial separa- 
tion of the holes D/M, where M is the Arnowitt-Deser- 
Misner (ADM) mass of the system. 

As an alternative to these two conformally fiat data 
types, the Lean code allows the iiac of nonspinning black- 
hole binary data based on the Kerr-Schild solution for a 
single black hole [68, 69]. The invariance of the structure 
of the Kerr-Schild data under boost transformations has 
motivated their use in boosted, superposed form. Even 
though these superposed data do not exactly satisfy the 
Einstein constraints for finite separation of the holes, 
they have been studied extensively in the literature, both 
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as initial data and in the context of binary-black-hole 
evolutions (see, for example, [51, 70-74]). Whenever we 
speak of superposed Kerr-Schild data in the remainder 
of this work, we thus refer to this direct superposition of 
the data which has been used in evolutions before. We 
are currently not aware of any time evolutions presented 
in the literature that start with Kerr-Schild data after 
applying a constraint solving procedure. 

In order to construct the data, we follow the approach 
of Ref. [51]. The solution of a single boosted, nonspinning 
Kerr-Schild hole with mass parameter m and velocity 
is calculated according to the prescription presented in 
App. B. The two solutions thus obtained for black holes 
at positions ^a;* and ^a;' are then superposed according 
to 



KS 
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-1/2 



(7) 
(8) 
(9) 

(10) 



We note that with this specific superposition, lapse a 
and shift /?' obey the close-limit condition, i. e. they lead 
to the lapse and shift of a single Kerr-Schild hole in the 
limit of zero separation. 



C. Gauge conditions 

An important ingredient in the recent success of nu- 
merical simulations of black-hole binaries has been the 

implementation of improved gauge conditions. In terms 
of the "3-1-1" decomposition, the coordinate invariance of 
general relativity is represented by the freedom to arbi- 
trarily specify the lapse function a and the shift vector /3' . 
While the particular choice of these functions leaves un- 
affected the physical properties of the spacetime, it can 
have a dramatic effect on the stability properties of a 
numerical simulation. 

In the past, the majority of gauge conditions have been 
designed with the purpose to drive the system of variables 
towards a stationary configuration (see e. g. [16, 75, 76]). 
In combination with the use of comoving coordinates, 
this approach lead to the first simulation of a complete 
binary-black-hole orbit [17]. More recent developments, 
however, have shown a tendency towards allowing the 
black holes to move across the computational domain 
(see e. g. [51, 53] for single moving black holes and head- 
on collisions and [15, 18, 19, 23] for orbiting black holes). 
We have implemented in the Lean code both the use 
of algebraic gauge conditions along the lines reported 
in [53] as well as live-gauge conditions similar to those 
presented in [18, 19] for the evolutions of black holes of 
the moving-puncture type (see also [77] for a more de- 
tailed numerical study and [78] for an analytic study of 
these types of gauge choices). Experimentally, we have 
found variations in these live-gauge conditions to mani- 
fest themselves most conspicuously in the profile of the 



variables near the punctures. In particular, we have 
noticed that steep gradients in these functions resulted 
in poor convergence properties of the merger time of the 
black holes, or, worse, instabilities. We have found opti- 
mal performance of our code in this respect by evolving 
the gauge variables according to 



dta = (3'dia-2aK, 



(11) 
(12) 
(13) 



Initially we have experimented with setting rj = 2, 
but observed an instability in the outermost refinement 

boundary for coarser resolutions. We have found the 
choice 77 = 1 to cure that instability while preserving the 
good convergence properties of the code and therefore use 
this value throughout this work. The gauge variables are 
initialized by using zero shift with a precoUapsed lapse 
a = = 

The gauge conditions (11)-(13) do not only provide 
stable evolutions, but also facilitate a comparatively sim- 
ple method to track the black-hole position. As has been 
shown in Ref. [24] the vanishing of x at the puncture in 
conjunction with Eq. (A6) implies that 



(14) 



We have implemented this relation via interpolation of 
the shift vector at the puncture location and subsequent 
update of the position using a second-order Runge-Kutta 
method. In practice we find excellent agreement between 
the resulting locations of the puncture and the coordinate 
center of the apparent horizon as calculated by AHFind- 
erDirect from surface integrals of the global coordi- 
nates over the horizon. 

In the case of the evolutions of Kerr-Schild data, we 
have also experimented with these gauge conditions. So 
far, however, we have not managed to obtain long-term 
stable evolutions in this way. We have therefore reverted 
to the approach of using algebraic gauge according to the 
procedure described in [51]. That is, we prescribe ana- 
lytic trajectories ^a;'(i), ^x^{t) for black holes A and B 
and calculate the resulting gauge functions by superpos- 
ing the analytic gauge of the individual holes. Following 
[53] we prescribe the analytic slicing condition in the form 
of the densitized lapse Q. We thus obtain 



KS^i = KS^^.^(A^^.+B^.)^ 



KS/ 



2 + 



1) 



-1/2 



(15) 
(16) 



Here the quantities denoted with an A or B are the ana- 
lytic expressions for the individual black holes and ^^7ij 
is the superposed metric defined in Eq. (7). In practice, 
we calculate the lapse from its densitized counterpart 

and the determinant of the numerical three-metric via 
1/2 

a = TnumQ. We emphasize that we use the densitized 
lapse only for algebraic slicing, but work with the un- 
modified lapse in all simulations using live-gauge condi- 
tions. 
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The trajectories used to evaluate the positions and ve- 
locities for the gauge functions associated with the in- 
dividual black holes are obtained from fifth order poly- 
nomials + vH + + ft^l<o -t- qH'^I^A. during the 
earlier stages of the infall of the black holes. In a time 
interval t\ < t < we perform a smooth (up to the 
fourth derivative) transition of these polynomials to the 
static function x* = 0. By virtue of the close-limit prop- 
erty of the superposed gauge (15), (16), we thus obtain 
a smooth transition of the gauge to that of a single non- 
spinning Kerr-Schild hole. 

The most difficult part in this procedure is to deter- 
mine the coefficients , a*, j% and fi, ti so that one 
obtains a stable simulation. We have only managed to 
obtain stable evolutions using gauge trajectories that are 
close to the coordinate trajectories of the apparent hori- 
zon center, in particular in the late stages of the infall. In 
practice the black holes collide along the z axis and we set 
the X and y components of all coefficients to zero. The re- 
maining z components are then obtained itcratively: The 
black holes are evolved with some initial guess for the co- 
efficients (normally those used for the simulation with the 
next smaller initial separation). The apparent horizon 
center is tracked until this simulation becomes unstable 
and we adjust the parameters to make the gauge trajec- 
tory agree better with the horizon motion. This process 
is repeated until a stable simulation is obtained. Some 
minor variations of the parameters are possible while pre- 
serving the stability of the simulation but do not have a 
significant impact on the resulting waveforms as is dis- 
cussed below in Sec. IV B 2. The exact parameters used 
for the Kerr-Schild simulations in this work are given be- 
low in Table III. Unless specified otherwise, we use the 
trajectories labelled "a" in that table. 



D. Black-hole excision 

Evolutions of puncture- type initial data have been per- 
formed in the past both with and without the use of 
black-hole excision (see e.g. [16, 17, 21, 48, 75]). Those 
without excision have commonly been achieved by factor- 
ing out the irregular part of the conformal factor while 
evolving only the regular remainder. It is a remarkable 
and surprising feature of the moving-puncture evolutions 
introduced in [18, 19], that these evolutions have been 
successful using neither excision, nor the factoring out 
of the irregular part of the conformal factor. Below we 
will follow the same approach for our puncture/Brill- 
Lindquist and Misner evolutions. 

In order to evolve Kerr-Schild data, however, we need 
to use black-hole excision. In contrast to puncture data, 
the spatial slices of the Kerr-Schild data do contain the 
physical singularity of the black hole at r = 0, which 
needs to be removed from the computational domain. 
Inside the Lean code we have implemented black-hole 
excision using either one-sided derivatives or extrapola- 
tion techniques. So far, we have obtained better stabil- 



ity properties using extrapolation which is the method 
of choice for all simulations presented in this work. This 
particular excision algorithm has been described in detail 
in Refs. [51-53]. In the Lean code the moving excision 
has now been generalized to work with moving refine- 
ment components. For this purpose each black hole has 
been assigned a particular refinement level it resides in 
(the finest level in all simulations presented in this work) . 
Excision for this black hole is then only performed on this 
refinement level and communication to coarser levels is 
performed exclusively via the restriction procedure inside 
Carpet. Special care must be taken in the black-hole 
excision if the refinement component has been moved be- 
cause the integer grid indices i, J, A; no longer correspond 
to the same coordinate position x,y,z as on previous 
time steps. Because the set of excision boundary points 
is stored in terms of their indices i,j, k, rather than their 
coordinate positions, we must recalculate the list of exci- 
sion boundary points every time the refinement compo- 
nent moves. This process docs not involve changing any 
of the BSSN variables, however; it merely corrects the 
book-keeping of the excision mask. 

With a correct list of excision points available at every 
time step we thus apply extrapolation of the BSSN vari- 
ables via second-order polynomials during each iteration 
of the Iterated Crank- Nicholson (ICN) cycle according to 
the procedure in Sec. 3 of Rcf. [52]. After the completion 
of the whole time step, the code checks for the position of 
the black hole and adjusts the center of the excision re- 
gion if necessary. As a minor modification compared with 
the excision method of the Maya code used in [51 53], 
we use the horizon finder to track the black-hole motion 
and move the excision region accordingly. 

So far, we have not succeeded in combining black-hole 
excision with the fourth-order discretization of the spa- 
tial derivatives. The problems largely arise from the need 
to use an excision boundary of thickness > 2 to accom- 
modate the wider fourth-order accurate stencils. For this 
reason, we use second-order discretization in space for all 
simulations using black-hole excision. 



E. Mesh refinement 

A further area of remarkable progress in numerical rel- 
ativity in recent years is that of mesh refinement, which 
is used almost routinely now in various forms in black- 
hole simulations. The need for using mesh refinement 
or essentially equivalent techniques based on specially 
adapted coordinates such as the "fish-eye transforma- 
tion" [79], arises from the presence of vastly different 
length scales in the spacetimes. On the one hand, a 
code has to resolve the steep gradients near the black- 
hole horizon, typically leading to length scales compara- 
ble with the mass of the hole. On the other hand, the 
typical wavelength associated with the ringing of a black 
hole is one order of magnitude larger. Furthermore, the 
calculation of accurate waveforms makes it necessary to 
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extract waves at sufficiently large radii, ideally, in the 
wave zone. This requires the use of computational grids 
at least two orders of magnitude larger than the radius of 
a single black hole. With current computational power, 
this can only be achieved inside the framework of mesh 
refinement. Simulations of moving black holes add the 
extra requirement of dynamic or adaptive refinement. 

In the Lean code, mesh refinement is provided by the 
Carpet package. Dynamic refinement based on Car- 
pet has already been reported in [80] . Here we use a re- 
fined version of this method. Carpet provides a routine 
which performs a regridding operation at regular inter- 
vals. That is, it interpretes a steerable parameter string 
which contains the exact specifications of all refinement 
components in terms of their corner positions. Inside the 
Lean code, we control this parameter string via a sep- 
arate thorn RegridInfo which works as follows. This 
thorn creates a map between each refinement component 
and the black hole it is tied to (a zero entry meaning 
the component is not tied to black-hole motion and re- 
mains stationary). The black- hole motion, in turn, is 
monitored, either using the horizon finder or the punc- 
ture tracking method according to Eq. (14). The cor- 
ner positions of the refinement components are adjusted 
according to the motion of the black holes. The Re- 
GRIdInfo thorn further performs checks on the internal 
consistency of the grid specifications and, if necessary, 
expands a component to guarantee that all finer compo- 
nents are accommodated with a minimum number of grid 
points between the refinement boundaries. Similarly, it 
expands components once the black-hole position comes 
too close to a refinement boundary. 

Finally, the thorn allows for the merger of previously 
separate refinement components. This is triggered by 
the distance between two components decreasing below 
a user-specified threshold value. Again, the parame- 
ter string used by Carpet is updated accordingly and 
the regridding completes the dynamic adjustment of the 
mesh refinement. We find this technique to work very re- 
liably and to preserve remarkably well the expected con- 
vergence properties of the code, as will be demonstrated 
below in Sec. III. 

For a given simulation the initial grid consists of two 
types of cubic refinement levels, n outer levels centered 
on the origin which remain stationary throughout the 
simulation and m levels with two components centered 
around either black hole. In the remainder of this work 
we specify the exact setup by giving the resolution h on 
the finest level as well as the radius of the cubes exclud- 
ing ghost zones required for interprocessor communica- 
tion. The grid spacing always increases by a factor of two 
from one level to the next coarser refinement level. For 
example, 

{(256, 128, 74, 24, 12, 6) x (1.5, 0.75), h = 1/48} 

specifies a grid with six fixed outer components of radius 

256, 128, 74, 24, 12 and 6 respectively and two refinement 
levels with two components each with radius 1.5 and 0.75 



centered aroimd either hole. The resolution is h = 1/48 
on the finest level and successively increases to 8/3 on 
the outermost level. In this work we will use equatorial 
as well as octant symmetry which reduces the number of 
points by a factor of 2 or 8 respectively. The grid setups 
used for the simulations of this work are summarized in 
Table I. 



F. Discretization of the BSSN equations 

In App. A we have listed explicitly the (f> and x version 
of the BSSN equations as used in the Lean code. The 
discretization of the spatial derivatives has been imple- 
mented in the form of second-order as well as fourth-order 
accurate stencils. With the exception of the advection 
derivatives of the form /3*(9i/ these stencils are centered. 
Advection derivatives, on the other hand, are approxi- 
mated with lop-sided stencils 

dxf = {—fi+2di,j,k + 4:fi+di,j,k — ^fi,j,k) , (17) 

dxf = Y2dx ^•^*+'^''*'-''*' ~ ^fi+2di,j,k + ^S>fi+di,j,k 

~10/ij,fc — 3fi-di,j,k) 7 (18) 

respectively for second and fourth-order accurate dis- 
cretization, with di = sgn{l3^) and likewise for the y and 
z direction. 

Using these representations for the spatial derivatives, 
the partial differential equations for the BSSN variables 
are integrated in time using the method of lines. Here the 
time discretization is performed using either the second- 
order accurate iterated Crank-Nicholson (ICN) scheme 
with two iterations [81] or standard fourth-order Runge- 
Kutta (RK) integration. The exact numerical implemen- 
tation of the BSSN equations is thus determined by three 
parameters, the time integration scheme, the ^ or x ver- 
sion and the order of the spatial discretization. In the re- 
mainder of this work, these different choices arc referred 
to as RKx4, ICN02 and so on. The discretizations used 
for the simulations in this work are summarized together 
with the grid setups in Table I. 

The Berger-Oliger-type mesh refinement provided by 
Carpet requires communication between the refinement 
levels, the so-called prolongation and restriction opera- 
tion (see e. g. [55] ) . For fourth (second)-order discretiza- 
tion of the spatial derivatives we use sixth (fourth)-order 
accurate prolongation in space with a total of nine or 
three buffer zones respectively for the RK and ICN time 
discretization (cf. Sec. 2.3 of [55]). Prolongation in time is 
always performed using third order accuracy. The neces- 
sary infrastructure required for higher order prolongation 
in time and, thus, genuine fourth-order accurate commu- 
nication between the refinement levels is not available 
in the currently used implementation of the mesh re- 
finement. The fourth-order convergence found for the 
simulations of puncture and Brill-Lindquist data below 
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TABLE I: Grid setup and numerical schemes used for the 
simulations presented in this work. The resolutions used for 
the convergence studies are hi = 1/48, /i2 = 1/44, ^3 = 1/40 
for models Rl and BL2, hi = 1/28, hi = 1/24, /ig = 1/20 
for model KS4 and hi = 1/400, /12 = 1/360, hs = 1/320 for 



model M4. 


Model 


Scheme grid 


Rl 


ICNx4 {(192, 128, 74, 24, 12, 6) x (1.5, 0.75), hi} 


BLl 
BL2 
BL3 
BL4 

ISCO 


RKx4 {(256, 128, 96, 32, 16) x (4, 2, 1), 1/48} 
RKx4 {(256, 128, 96, 32, 16) x (4, 2, 1), hi} 
RKx4 {(256, 128, 96, 32, 16) x (4, 2, 1), 1/48} 
RKx4 {(256, 128, 96, 32, 16) x (4, 2, 1), 1/48} 
ICN04 {(256, 128, 88, 24, 12, 8) x (2.4, 1.2, 0.6), 1/40} 


KSl 
KS2 
KS3 
KS4 


ICN</>2 {256, 128, 96, 32, 16) x (4, 2), 1/24} 
ICN</>2 {256, 128, 96, 32, 16) x (4, 2), 1/24} 
ICN<^2 {256, 128, 96, 32, 16) x (4, 2), 1/24} 
ICN<^2 {256, 128, 96, 32, 16) x (4, 2), hi} 


Ml 
M2 
M3 
M4 


ICNX4 {60,30,f,f,f)x(l|,i|,l|),^} 
ICNX4 {48,24,18,6,3) X (|,|,^),_3_ 
ICNX4 {40,20,15,5,1) X (f, A, 
ICNX4 {32, 16, 12.8, 4, 2) x (i, i, i), h,} 



indicates that this does not represent a problem for the 
type of simulations under discussion in this work. 



G. Wave extraction 

We extract gravitational waves from our numerical 
simulations by calculating the Newman-Penrose scalar 
^'4 using the electromagnetic decomposition of the Weyl 
tensor which is described in Appendix C. The spatial 
derivatives required in this calculation are obtained us- 
ing either second or fourth-order accurate stencils cho- 
sen in accordance with the spatial discretization of the 
BSSN evolution equations. The calculation of ^4 as well 
as the extraction of modes has been tested with the an- 
alytic expression calculated for the Teukolsky wave [82] 
for both the £ = 2, m — and £ — 2, m = 2 wave [cf. also 
Ref. [83]]. In both cases, the evolutions have been car- 
ried out using the ICN(/)2 implementation of the BSSN 
equations and resulted in second-order convergence of the 
waveforms. 

Once ^'4 has been calculated on a sphere of constant 
extraction radius, the radiated energy and momenta are 
obtained from Eqs. (22)-(24) in Ref. [84]. In practice we 
perform these calculations in a post-processing operation 
using the output data of '^4. There, we calculate both 
the total radiated energy as well as the energy radiated 
in the dominant modes, £ — 2, m ~ ±2 for orbiting con- 
figurations and £ = 2, m = for the head-on collisions. 
For all head-on collisions we find the dominant mode to 
be responsible for > 99% of the total radiated energy; 
for the inspiral the dominant modes account for about 



98.5 % of the total energy. 

III. BINARY BLACK-HOLE ORBITS 

Before we compare the head-on collisions of differ- 
ent data types, we demonstrate the code's capability to 
produce evolutions of orbiting black-hole binaries with 
convergent waveforms. For this purpose we consider 
model Rl of Table I of Ref. [19]. Here two black holes 
with mass parameter m = 0.483 start at coordinate po- 
sitions X = ±3.257 with linear momentum parameter 
P = ±0.133 in the y direction. 

We evolve this configuration with a setup as specified 
for model Rl in Table I using three different resolutions 
hi = 1/48, h2 = 1/44 and /13 = 1/40. The simula- 
tions are performed using equatorial symmetry across the 
orbital xy plane. 

The resulting real part of the £ = 2, m — 2 mode of 
the Newman-Penrose scalar ^'4 extracted at r = 60M 
is shown in Fig. 1 for all three resolutions. We first note 
that the waveforms show good agreement with the results 
obtained from similar simulations in the literature [27, 
35]. A factor two discrepancy with Fig. 2 of [27] results 
from a trivial rescaling depending on the choice of the 
eigenmode basis [cf. their Eq. (4)]. 

With regard to a convergence analysis, we first note 
that the error manifests itself in two forms, a phase shift 
and an amplitude difference. We therefore study the con- 
vergence both with and without applying a phase cor- 
rection to align the global maxima of the curves. For 
the convergence analysis one commonly assumes that 
the discretization error be dominated by a leading or- 
der term 0{h°'), so that the numerical solution fh of a 
given grid function is related to the continuum limit / by 
fh ~ f + Ch", where the coefficient C does not depend 
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FIG. 1: Real part of the ^ = 2, m = 2 multipole of Mr*4 
extracted from the Rl simulation at rex ~ 60Af obtained for 
resolutions hi, /12 and hs. 
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FIG. 2: Convergence analysis of the £ — 2, m = 2 multipole of 
without correcting the phase error (upper panel) and 
after applying a phase correction (lower panel). 



50M, 60M and 70 M. We can use these results to esti- 
mate the uncertainties in the radiated energy resulting 
from finite resolution and extraction radii. The stan- 
dard procedure to assess the impact of the resolution is 
to apply Richardson extrapolation, i. e. extrapolate the 
values obtained for a convergent simulation to the con- 
tinuum limit /i — > 0. Using this procedure we obtain 
Etot = 3.558 %, 3.543 % and 3.532 % respectively of 
the total ADM mass M of the system at extraction radii 
50M, 60M and TOM. This corresponds to an estimate of 
the discretization error in the radiated energy of about 
1 % for the high resolution hi = 1/48 and about 2 % for 
the low resolution hs = 1/40. 

In complete analogy to the procedure used to study the 
convergence with grid resolution h, we use these values 
to estimate the dependency of the radiated energy on 
the extraction radius. We find the resulting error to be 
modeled well by a l/vcx fall-off, i. e. 

^^tot « Etotl^^^^ + O (^-^^ . (20) 



m 




10 

(x+l)/M 



FIG. 3: Convergence analysis of the Hamiltonian constraint 
on the X axis at t = 128 M, shortly after the crossing of the 
punctures. 

on the resolution h. Applying this relation to the three 
different resolutions /ii, /i2, ^3 one obtains 

fhj - fh2 _ ^3 ~ ^2 j--|^g-j 

fh2 - fhi - /l" 

Applied to our case, this relation leads to the value 
1.58 for the case of fourth-order convergence. The con- 
vergence behavior of the £ = 2, to 2 multipole of 
with and without a phase correction is shown in 
Fig. 2, where we have amplified the differences between 
the higher resolution runs by the factor 1.58 expected 
for fourth-order convergence. The analysis shows good 
agreement with fourth-order convergence in both cases. 

We similarly observe fourth-order convergence for the 
total radiated energy extracted at coordinate radii = 



Extrapolation of the results obtained at finite extraction 
radii thus gives a total radiated energy of E'tot = 3.466 % 
of the total ADM mass as well as an estimate for the 
error arising out of the use of a finite extraction radius 
of 2.7 %, 2.3 % and 2.0 % respectively for r^^ = 50M, 
60M and TOM. 

For the simulations presented in this work we univer- 
sally find the errors due to finite differencing and finite 
extraction radius to point in opposite directions: finite 
resolution leads to underestimating the amount of radi- 
ated energy, a finite extraction radius overestimates the 
energy. We therefore feel justified in using the sum of the 
individual errors as a conservative upper limit for the to- 
tal error. In this case, we obtain a numerical error of 3 % 
for the high resolution simulation using = 70M . 

Repeating the same calculation without including the 
artificial radiation burst due to the initial data^, we ob- 
tain a total radiated energy of iJ^ad = 3.408 % of the 
ADM mass. We note that this result for the energy shows 
excellent agreement with those presented in the literature 
(cf. Table III of [27]). 

As a further test of the code we follow Rcf. [18] and 
check the convergence properties of the Hamiltonian con- 
straint on an equatorial axis shortly after the crossing 
of the punctures. The result is shown in Fig. 3 where 
we have amplified the high resolution result by a fac- 
tor of 1.2^ as expected for fourth-order convergence. In 
spite of the presence of some numerical noise, the fig- 
ure demonstrates compatibility with overall fourth-order 
convergence of the simulations. We believe the larger 
amount of noise, as compared with the results of [18], to 



^ In practice we ignore contributions at t < rex + 30A/ in the 
waveforms. 
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be a consequence of the use of mesh refinement and the 
discontinuous error terms at the refinement boundaries. 

Unfortunately we are currently not able to obtain sim- 
ilar orbital simulations with Kerr-Schild data for want 
of suitable live-gauge conditions analogous to Eqs. (11)- 
(13). We therefore perform the comparison between these 
two data types and the Misner data inside the framework 
of head-on collisions. 



IV. HEAD-ON COLLISIONS 

Head-on collisions represent the simplest form of black- 
hole binaries and have been studied numerically in var- 
ious forms for a long time. The majority of such sim- 
ulations has been performed using data of Misner [67] 
or Brill-Lindquist type (see e.g. [16, 83, 85-88]). As an 
alternative, collisions using Kerr-Schild data have been 
investigated in [51]. Here we will study in detail head-on 
collisions of all three data types and compare the results. 

The time evolutions of these two types of initial data 
present different difficulties and therefore require differ- 
ent evolution techniques. In summary, these are the use 
of second-order differencing, algebraic gauge conditions 
and black-hole excision for the Kerr-Schild data, whereas 
Brill-Lindquist data are evolved using fourth-order dis- 
cretization without using the black-hole excision proce- 
dure described in Sec. HD. As will be demonstrated 
below, we have obtained satisfactory convergence perfor- 
mance for the Kerr-Schild evolutions using the </> version 
of the BSSN equations. In the case of Brill-Lindquist 
data, we find the x version more successful in providing 
fourth-order convergence. The Misner data are conceptu- 
ally similar to Brill-Lindquist data. Experimentally, how- 
ever, we have found the Misner data to lead to substan- 
tially larger amounts of numerical noise originating from 
the refinement boundaries when evolved in time with the 
Runge-Kutta method. Below we will show that the noise 
level is acceptable when using the second-order accurate 
ICN scheme instead. 

In the notation of Sec. II F we therefore evolve the 
Kerr-Schild data using the ICN02 scheme, the Brill- 
Lindquist data with the^ RKX4 and Misner data with 
the ICNx4 scheme. The resulting accuracy and conver- 
gence properties will be studied in detail in Sec. IV B. 



A. Choice of initial parameters 

A fundamental difficulty in the comparison between 
simulations of Brill-Lindquist, Kerr-Schild and Misner 
data is the physical interpretation of the initial-data sets. 



^ with the exception of the simulations in Figs 7 and 13 which are 
not used in this quantitative comparison 
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FIG. 4: Binding energy Et/M for Kerr-Schild (solid), Brill- 
Lindquist (short-dashed) and Misner (long-dashed) initial- 
data sets as function of the coordinate distance D/M of the 
holes. 

We first note that there exists no general method to rig- 
orously quantify the degree to which two such initial con- 
figurations represent the same physical scenario. As an 
approximation, we determine the initial parameters as 
follows. First, we start the head-on collisions with two 
black holes of equal mass at rest, and thus eliminate the 
question of choosing initial linear momenta and mass ra- 
tios. 

Except for a rescaling of the entire spacetime corre- 
sponding to a rescaling of the system's total mass, an 
initial configuration for a head-on collision of nonspin- 
ning, equal-mass black holes is characterized by one free 
parameter which can be viewed as a measure for the ini- 
tial separation of the black holes or the binding energy 
of the system. In the case of Misner data, this degree 
of freedom is represented by the parameter /i in Eq. (5), 
whereas Brill-Lindquist and Kerr-Schild data require co- 
ordinate positions of the black holes so that the free pa- 
rameter is the initial coordinate separation D. For the 
comparison of the different data types we fix the remain- 
ing free parameter by demanding that all three versions 
of the binary-black- hole spacetime have identical binding 
energy 

M ' ^ > 

where the irreducible black- hole masses Mi , M2 are given 
by their respective apparent horizon masses. For illustra- 
tion we plot in Fig. 4 the binding energy as a function of 
the coordinate distance D. For Misner data, the initial 
black-hole centers are approximated by their apparent 
horizon position which is ±1.0 for all simulations dis- 
cussed in this work. 

The freedom in rescaling the spacetimes is fixed by 
demanding the total ADM mass to be unity for Kerr- 
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Schild and Brill-Lindquist data which implies bare mass 
parameters nii — m2 — 0.5 in both cases. In contrast, 
the conformal factor for the Misner data in Eq. 5 does 
not contain a bare mass parameter and the ADM mass 
depends on the value of /i. Specifically, it decreases for 
larger fi and thus implies a larger black-hole separation 
D/M. In our simulations we have taken into account the 
different ADM masses of Kerr-Schild and Brill-Lindquist 
data on the one hand and Misner data on the other by 
using numerically higher resolutions for the Misner mod- 
els. Relative to the ADM mass, however, the resolutions 
are rather similar to those used for Brill-Lindquist data. 
The exact ADM masses for all simulations discussed in 
this comparison are given in the third column of Table 
II. 

Using this procedure, we have determined four dif- 
ferent models with binding energies Eh/M — —0.0290, 
-0.0240, -0.0197 and -0.0169. These models are listed 
in Table II as BLl - 4, KSl - 4 and Ml - 4. 



B. Testing the code 

In this section we calibrate the performance of the code 
in the case of head-on collisions of all data types by per- 
forming convergence tests and investigating other error 
sources. For Brill-Lindquist data we also use results pub- 
lished in Ref. [83] starting from the approximate separa- 
tion of the innermost stable circular orbit (ISCO). 

1. Brill-Lindquist data 

We assess the discretization error of the Brill-Lindquist 
simulations by evolving model BL2 of Table II using three 
different resolutions hi = 1/48, /12 = 1/44 and /13 = 1/40 
with a constant Courant factor of dt/dx = 1/2. We 
have studied the convergence of the resulting gravita- 
tional waveforms in complete analogy to the procedure 
used in Sec. Ill for black-hole inspiral. We observe fourth- 
order convergence as is demonstrated in Fig. 5 for the 
£ = 2, m = multipole extracted at = 40M. We 
similarly observe fourth-order convergence for the total 
radiated energy and use Richardson extrapolation to es- 
timate the discretization error. We find the relative error 
at a resolution h — 1/44 to be less than 0.5 % and use 
this value as a conservative upper limit. 

In order to assess the impact of extracting waves at 
finite radii we have studied the wave signal at extraction 
radii 40M, TOM and 90M. The resulting £ = 2, m = 
multipole of the Newman-Penrose scalar ^'4 is shown in 
the upper panel of Fig. 6. We estimate the uncertainty 
due to the extraction radius in the same way as in Sec. Ill 
and find the relative error in the total radiated energy to 
be of the order of 1 % for rex = 40M and less for the 
larger radii TOM and 90M. 

As a further test of our code, we compare the wave- 
forms obtained for Brill-Lindquist data with those avail- 
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FIG. 5: Convergence analysis for the I = 2, m = Q mode of 
the Newman-Penrose scalar '^4, extracted at rex = 40M. The 
difference between the runs with higher resolutions has been 
amplified by a factor 1.58 expected for fourth-order conver- 
gence. 
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FIG. 6: The £ = 2, m = multipole obtained for models BL2 
(upper) and KS2 (lower panel) extracted at rex = 40M, TOM 
and 90M. 

able in the literature. The head-on collisions presented 
commonly start with time symmetric initial data of two 
holes at positions ±1.1515 Af. This value has been cal- 
culated in Refs. [89, 90] for the ISCO. The £ = 2, to = 
mode of the Newman-Penrose scalar ^'4 of this configu- 
ration has been calculated in [83] at an extraction radius 
Tex = 20M. In Fig. T we plot our result extracted at 
the same radius and obtained using the setup labelled 
"ISCO" in Table I. Up to the trivial rescaling factor of 
two mentioned above, we find excellent agreement with 
Fig. 5 in Ref. [83]. To our knowledge, similar results ob- 
tained with larger initial separations of the holes have not 
yet been published. It is part of the motivation of this 
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TABLE II: Summary of the simulations performed in this work. Simulation Rl is the inspiral simulation described in Sec. III. 
The other simulations are the head-on collisions performed for the comparison of Brill-Lindquist, Misner and superposed Kerr- 
Schild data. M is the total ADM mass of the spacetime, D the initial coordinate separation of the holes (for Misner data we list 

the parameter /j, instead) and £5 the binding energy M — Mi — A/2- -Etot, -Eini and E^ad are the total radiated energy, the energy 
contained in the spurious initial burst and the energy radiated in the inspiral and merger. The uncertainties included are those 
due to finite differencing, finite extraction radius and the uncertainties in separating the merger signal from the spurious iuitial 
burst. For the Kerr-Schild data we also list the radiated energies obtained from extrapolation to rex — > 00 with uncertainties 
due to finite differencing and the interference of the initial burst. 



Model 


D ov n 


M 




Etot/M 


Eini/M 


Erad/M 


Rl 


6.5 


0.996 


-0.0145 


(3.466 ±0.104) % 


(0.058 ± 0.002) % 


(3.408 ±0.102) % 


BLl 


8.6 


1 


— n 02Q0 


(0 0553 + n 0008"! % 


(0 0031 + 0001 "1 % 


(n 0522 + 0008") % 


BL2 


10.2 


1 


-0.0240 


(0.0553 ± 0.0008) % 


(0.0022 ± 0.0001) % 


(0.0531 ± 0.0008) % 


BL3 


12.5 


1 


-0.0197 


(0.0557 ± 0.0008) % 


(0.0014 ± 0.0001) % 


(0.0543 ± 0.0008) % 


BL4 


14.6 


1 


-0.0169 


(0.0564 ± 0.0009) % 


(0.0009 ± 0.0001) % 


(0.0555 ± 0.0008) % 


KSl 


10.0 


1 


-0.0290 


(0.1099 ±0.0175) % 


(0.0540 ± 0.0119) % 


(0.0560 ± 0.0123) % 


Tex — > 00 








(0.0963 ± 0.0029) % 


(0.0438 ± 0.0049) % 


(0.0525 ± 0.0052) % 


KS2 


12.0 


1 


-0.0240 


(0.0962 ± 0.0154) % 


(0.0325 ± 0.0072) % 


(0.0617 ±0.0119) % 


rex 00 








(0.0844 ± 0.0025) % 


(0.0284 ± 0.0019) % 


(0.0560 ± 0.0032) % 


KS3 


14.0 


1 


-0.0197 


(0.0888 ± 0.0142) % 


(0.0227 ± 0.0040) % 


(0.0661 ± 0.0110) % 


rex — » 00 








(0.0789 ± 0.0024) % 


(0.0198 ± 0.0010) % 


(0.0591 ± 0.0023) % 


KS4 


16.0 


1 


-0.0169 


(0.0855 ± 0.0137) % 


(0.0163 ± 0.0028) % 


(0.0692 ± 0.0112) % 


rex 00 








(0.0751 ± 0.0023) % 


(0.0140 ± 0.0007) % 


(0.0611 ± 0.0021) % 


Ml 


3.573 


0.231 


-0.0290 


(0.0555 ± 0.0017) % 


(0.0031 ± 0.0001) % 


(0.0524 ± (0.0016) % 


M2 


3.757 


0.191 


-0.0240 


(0.0556 ± 0.0017) % 


(0.0021 ± 0.0001) % 


(0.0535 ± (0.0016) % 


M3 


3.948 


0.157 


-0.0197 


(0.0560 ± 0.0017) % 


(0.0014 ± 0.0001) % 


(0.0546 ± (0.0017) % 


M4 


4.096 


0.135 


-0.0169 


(0.0567 ± 0.0017) % 


(0.0009 ± 0.0001) % 


(0.0558 ± (0.0017) % 
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FIG. 7: The £ = 2, m = mode of Mr*4 extracted at 
rex = 20M from a head-on collision of Brill-Lindquist data 
starting from the approximate ISCO separation (cf. Fig. 5 in 
Ref. [83]. 



work to provide such an extension of the existing work. 
We finally address one conceptual difference between 

the Brill-Lindquist and Misner data on the one side and 
Kerr-Schild data on the other. Whereas the former ini- 



tial data are time-symmetric, the superposed Kerr-Schild 
data do not satisfy this requirement, even for vanishing 
velocity parameter. Wc thus cannot rule out that the in- 
dividual Kerr-Schild holes do actually have a small boost 
and thus represent a slightly different physical configu- 
ration. Unfortunately, there exists no rigorous way to 
quantify the linear momenta of the individual holes in 
the Kerr-Schild spacetime, although the hypothesis of 
small momenta is compatible with the small nonzero ini- 
tial coordinate velocities of the apparent horizon posi- 
tions of the order of v = 0.05 towards each other which 
we observe for the Kerr-Schild holes. For this reason we 
proceed differently and instead consider the impact of 
small initial boosts as an additional uncertainty in our 
study. We quantitatively study this effect using a mod- 
ified version of model BL4. Specifically, we use punc- 
ture data, where initial linear momenta pointing towards 
each other arc applied to the individual black holes in the 
form of nonzero Bowen-York [66] momentum parameters 
P = 0.035 and P = 0.067. All other parameters for these 
puncture models are kept at the values of model BL4. 

In Fig. 8 we show the resulting waveforms at rex = 
70Af shifted in time to align their global extrema. The 
differences in the waveforms are rather small and we 
obtain for the energy radiated in the infall and merger 
E,^d = 0.05601±0.0008 and 0.05795±0.0009 respectively 
of the ADM mass. Compared with the nonboosted result 



12 



Mry\f. 




t/M 



FIG. 8: The I = 2, m = Q mode of Mr*4 extracted at 
»"cx = 70 Af for model BL4 (solid), as well as puncture initial 
data with P = =f0.0335 (dashed) and P = =f0.067 (dotted) 
for the holes starting at coordinate positions z — ±7.3. 
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FIG. 9: Convergence analysis of the £ = 2, m — mode of 
'ifi obtained at rex = 90M for model KS4 using resolutions 
hi = 1/28, h2 = 1/24 and /is = 1/20. The differences between 
the high resolution results has been amplified by a factor of 
1.66 expected for second-order convergence. 



0.0555±0.0008, this corresponds to systematic deviations 
of about 1 % and 4 %. 



2. Kerr-Schtld data 

By evolving the Kerr-Schild models of Table II at dif- 
ferent resolutions, we find model KS4 with the largest 
initial separation to exhibit the largest uncertainties. We 
therefore focus in our convergence analysis on this model 
to obtain conservative upper limits for the discretization 
error. We have evolved this model using the finest reso- 
lutions hi = 1/28, /i2 = 1/24 and /13 = 1/20 and a con- 
stant Courant factor of 1/4. Compared with the Brill- 
Lindquist data we have to use these seemingly coarser 
resolutions on the finest level because the coordinate ra- 
dius of the apparent horizon is larger in the Kerr-Schild 
case and the use of excision prohibits the use of refine- 
ment components inside the apparent horizon. We em- 
phasize, however, that relative to the coordinate radius of 
the horizon, about rah = 2Mi^2 for Kerr-Schild data and 
Tah = Mi^2/2 for Brill-Lindquist data, our setup results 
in a finer resolution in the Kerr-Schild case. 

The resulting differences in the I — 2, m — Q multipole 
of the Newman-Penrose scalar are shown in Fig. 9 and 
demonstrate second-order convergence. Using Richard- 
son extrapolation as before, we obtain an error of 3 % 
in the radiated energy due to the discretization of the 
Einstein equations. 

In comparison with the Brill-Lindquist data, we ob- 
serve larger differences in the amplitude of the gravita- 
tional waves extracted at different radii. This is shown in 
the lower panel of Fig. 6 where we plot the € = 2, m = 
multipole obtained for model KS2 at = 40Af, 70A/ 
and 90M. Systematically investigating the total radiated 



energy for all Kerr-Schild simulations we find the results 
compatible with a l/r^-^ fall-off as in Sec. III. For exam- 
ple, the energies extracted from model KS4 at 40Af, 70M 
and 90Af using a resolution ft, = 1/24 are 0.0984, 0.884 
and 0.855 % of the total mass M respectively. Extrapola- 
tion to infinite radius results in 0.0752 % of the mass and 
a relative error of 13 % at extraction radius r^x = 90Af . 
We find the uncertainties due to the extraction radius to 
be very similar for all other Kerr-Schild simulations and 
to be essentially independent of the grid resolution. 

In view of this large uncertainty we will always present 
in the remainder of this work two values for the energies 
resulting from simulations of Kerr-Schild data. The nu- 
merical values obtained at the largest extraction radius 
Tex = 90M, together with uncertainties due to extrac- 
tion radius, discretization and interference of the initial 
pulse, and the value extrapolated to infinite radius with 
uncertainties due to discretization and the initial pulse 
are listed in Table II. 

A further difficulty in the case of the Kerr-Schild data 
arises from the relatively large amount of spurious ra- 
diation due to the initial data. This spurious radiation 
manifests itself as a pulse in the waveform starting at 
t w Tex. In the two panels of Fig. 6, for example, the ini- 
tial burst leads to local extrema in the waveforms near 
t K, 50M. The amplitude of the pulse is, however, sub- 
stantially larger for the Kerr-Schild (lower panel) than 
the Brill-Lindquist (upper panel) data. For the smaller 
separations used in our analysis, it becomes nontrivial 
to disentangle this pulse from the actual merger signal 
and it is not entirely clear, how much radiated energy is 
due to the black-hole merger and how much due to the 
spurious pulse. We attempt to bracket these uncertain- 
ties by using a variable threshold ttiucsh so that radia- 



13 



TABLE III: Parameters for the gauge trajectories used for the 
Kerr-Schild simulations. 



Model 


z/M 




a" M 






ti 
If 


i2 
a7 


KSl a 


±5 





TO.037 


±0.0038 





10 


35 


Kbi b 


±5 


TO.08 tO.0061 


TO.0002 





20 


40 


KS2 a 


±6 





TO.029 


±0.004 


TO.000278 


25 


50 


KS2 b 


±6 


TO.06 


TO.008 


±0.0004 


TO.00002 


20 


44 


KS3 a 


±7 





TO.022 


±0.0027 


=F0.000165 


25 


57 


KS3 b 


±7 


TO.04 


TO.007 


±0.0003 


TO.000012 


25 


57 


KS4 a 


±8 





TO.018 


±0.002 


±0.000104 


34.5 


84.7 


KS4 b 


±8 


TO.03 


TO.006 


±0.00027 T0.000012 


50 


70 



tion at t < tthrosh is considered part of the initial pulse 
and radiation at larger t part of the merger signal. At a 
given extraction radius we then vary tthresh in the range 
rex±30M to rex±40M. The two resulting energy contri- 
butions are given by their average values obtained over 
this interval plus or minus an error given by the upper 
and lower bounds. 

We finally discuss the impact of the gauge trajectories 
on the resulting waveforms. We have already mentioned 
in Sec. II C that the gauge trajectories need to closely 
resemble the motion of the center of the apparent hori- 
zon and are obtained iteratively by approximating the 
horizon trajectory. For this comparison we have con- 
structed gauge trajectories according to this procedure, 
first, by fixing to be zero^ and, second, by also ad- 
justing this parameter in the iterative procedure. The 
parameters for the trajectories are listed in Table III and 
in Fig. 10 we show as examples the resulting waveforms 
extracted at = 90M for simulations KS2a, b (upper 
panel) and KS4a, b (lower panel) obtained using a reso- 
lution h = 1/24. The resulting waveforms are practically 
indistinguishable and the differences in the energies for 
the initial pulse, the merger signal and the total wave- 
form are about 1 % and thus significantly smaller than 
the uncertainties due to the finite differencing and the 
extraction radius. The same applies to variations from 
2/3 to zero in the evolution parameter a in Eq. A5. 



3. Misner data 

Finally, we estimate the numerical error of the evolu- 
tions starting from Misner data. All evolutions of these 
data using the Runge-Kutta time integration have re- 
sulted in significant contaminations of the resulting wave- 
forms by numerical noise. A comprehensive analysis of 
the performance of different numerical schemes in evolu- 



* This velocity parameter is not to be confused with that used in 
the calculation of the initial data for 7ij and Kij which is always 
zero. 
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FIG. 10: The ^ = 2, m = multipole of Mr*4 extracted 
at rex = 90M for model KS2 (upper) and KS4 (lower panel) 
using gauge trajectories labeled in Table III as "a" (solid) and 
"b" (dashed curve). 
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FIG. 11: Convergence analysis of the ^ = 2, m = mode 
of 'I'4 obtained for model M4 using resolutions h\ = 1/400, 
h2 = 1/360 and /13 = 1/320. The differences between the 
high resolution results has been amplified by a factor of 1.40 
expected for second-order convergence. 

tions of Misner data and the underlying causes is beyond 
the scope of this work, but we will show here that suffi- 
ciently accurate simulations can be obtained by using the 
ICN scheme instead of Runge-Kutta and also choosing a 
small Courant factor of 1/8. 

The resulting differences in the i = 2, m = mode 
of Afr^4 obtained for model M4 are shown in Fig. 11. 
Compared with the Brill-Lindquist and Kerr-Schild evo- 
lutions, we observe larger amounts of high frequency 
noise in the early stages of the simulations due to the 
spurious initial radiation. Still, the overall behavior is 
compatible with the expected second-order convergence. 
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Using the same methods as before, we find the resulting 
uncertainty in the radiated energy due to finite differ- 
encing to be of the order of 1 %. With regard to the 
extraction radius, we find Misner data to behave similar 
to Brill-Lindquist data. The resulting uncertainty due to 
the use of finite radii is about 2 % at r-ox — 40M and 1 % 
at TOM or 90M. 



C. Results 

In order to compare the different initial-data types, we 
have evolved all models listed in Table II with the grid 
setups described in Table I. We summarize the results in 
Fig. 12 which shows the £ — 2, m = modes obtained at 
extraction radius r^^ = 90M and Table II where we list 
the initial parameters and the total radiated energy i^tot 
as well as the contributions of the spurious initial pulse 
Eini and the merger signal i?rad- The radiated energies 
have been extracted at r^x = 90M and for Kerr-Schild 
data we also give the extrapolated values for rex — *■ oo. 
The uncertainties are those obtained in the previous sec- 
tion. In the case of Kerr-Schild data starting from small 
separations, the error for Ei^i and i?rad is amplified sig- 
nificantly by the uncertainties in separating the initial 
pulse from the merger signal. 

In Fig. 12 we note the substantially larger amount of 
artificial radiation due to the initial data in the Kerr- 
Schild case (dashed curve). Second, we observe excellent 
agreement between the waveforms obtained from Mis- 
ner and Brill-Lindquist data (the solid and long-dashed 
curves are practically indistinguishable) and good qual- 
itative agreement with the Kerr-Schild results. There 
remains, however, a small systematic deviation to the ef- 
fect that the Kerr-Schild waves have a 5 — 10 % larger 
amplitude. We have already seen, however, that the fi- 
nite extraction radius results in an overestimation of the 
wave amplitudes, in particular in the Kerr-Schild case. 

In order to quantify this effect, we consider the radi- 
ated energies in Table II. We first note that within the 
uncertainties Brill-Lindquist and Misner data result in 
identical amounts of energy in the initial-data pulse, the 
merger waveform as well as the total radiation. In con- 
trast, the total radiated energy obtained from the Kerr- 
Schild data is significantly larger than that of Misner 
and Brill-Lindquist simulations. This excess energy is 
largely due to the spurious initial pulse, however, and for 
models KSl and KS2, the energy contained in the physi- 
cally important merger waveform agrees within the error 
bounds with its Misner and Brill-Lindquist counterparts. 
In particular, this is true for the values obtained from 
extrapolation to infinite extraction radius. The situation 
becomes more complicated, however, for the Kerr-Schild 
models starting from larger separations, in particular for 
model KS4. With the error estimates obtained in the pre- 
vious section, we obtain a lower limit of E^iad = 0.0580 % 
of the ADM mass which exceeds the upper limit of sim- 
ulation BL4 of E,^d = 0.0563 by about 3 %. While we 
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FIG. 12: The 1 = 2, m = G mode of Mr*4 at r^x = 90M for 
the Brill-Lindquist, Kerr-Schild and Misner versions of models 
1-4. For presentation purposes, the Kerr-Schild data have 
been shifted in time to align the maxima of the waveforms. 



have taken into account in the derivation of these results 
the errors arising from finite differencing and the wave 
extraction at finite radius, it is possible that systematic 
errors are responsible for the remaining discrepancy. We 
have seen in Sec. IV B 1 that small boosts give rise to ra- 
diated energies a few per cent larger than those obtained 
for initially time symmetric data. A further systematic 
error results from the constraint violations inherent to 
the superposed Kerr-Schild data. Unfortunately there 
exists, to our knowledge, no literature on time evolu- 
tions of the constraint-solved version of the Kerr-Schild 
data. Filling this gap is beyond the scope of this paper 
as it requires the addition of elliptic solvers, currently 
not available in the Lean code. It is therefore currently 
not possible to rigorously quantify the impact of the con- 
straint violations on the resulting waveforms. We note, 
however, that the amount of spurious initial gravitational 
radiation inherent to the superposed Kerr-Schild data is 
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FIG. 13: Waveforms obtained from Brill-Lindquist data start- 
ing at various initial separations. The data have been shifted 
in time to align the maxima of the waveforms. 



significantly larger than the discrepancies we observe. If 
this spurious initial radiation is a signature of the con- 
straint violations, it is certainly possible that the discrep- 
ancies observed here are due to the constraint violations 
of the Kerr-Schild data. 

A further interesting question is the dependence of the 
radiated energies on the initial black-hole separations. 
We have already noticed that the amount of spurious 
initial radiation decreases at larger separations as is ex- 
pected. Correspondingly, we observe the expected in- 
crease in the energy radiated in the infall and merger 
at larger initial separations. This increase is relatively 
weak, though, for the cases studied here, especially for 
Brill-Lindquist and Misner data. It would be desirable 
to probe a larger range of initial distances, in particular 
smaller separations, to study this behavior in more detail. 
Unfortunately, such an extension encounters difficulties 
at either end of the spectrum. 

In the case of Kerr-Schild data, separations smaller 
than that of KSl lead to a severe contamination of the 
actual signal by the spurious wave content and thus do 
not allow a physically meaningful interpretation. At the 
upper end, we are limited by the construction of suitable 
gauge trajectories. The prolonged infall puts stronger de- 
mands on the fine tuning of the gauge trajectories. So far, 
we have not managed to obtain stable evolutions starting 
from Kerr-Schild data with D > IQ M. 

In the case of Brill-Lindquist data, we do not encounter 
such difficulties with the gauge because of the universal- 
ity of the live-gauge conditions (11)-(13). Results start- 
ing from small separations, however, are subject to the 
difhculties due to the initial wave burst. This is illus- 
trated by evolving a set of Brill-Lindquist data starting 
with initial separations Dbl — 2.3 M, A.3 M and 6.3 M. 
The resulting waveforms as obtained with the setup la- 
belled "ISCO" in Table I are shown together with that 



of model BLl in Fig. 13. We clearly notice a substantial 
contamination of the waveforms at small separations by 
radiation inherent to the initial data. In consequence, 
an accurate calculation of the gravitational wave energy 
generated in the head-on collision becomes highly non- 
trivial even for Brill-Lindquist data. For this reason the 
comparison performed in this study is currently limited 
to the window of binding energies or separations covered 
by the results in Table II. 



V. SUMMARY AND CONCLUSIONS 

In this work, we have presented in detail a numerical 
code designed for the simulation of black-hole binaries 
in the framework of three-dimensional, nonlinear general 
relativity. The code facilitates black-hole evolutions us- 
ing different initial-data types and evolution techniques. 

It has been demonstrated that the code is capable of 
evolving state-of-the-art binary-black-hole orbits using 
the recently developed moving-puncture technique. With 
regard to the accuracy of the results, we find it crucial 
to use a fourth-order discretization of the spatial deriva- 
tives appearing in the BSSN formulation of the Einstein 
field equations. The resulting simulations yield conver- 
gent waveforms which agree well with results presented 
in the literature. The same holds for the radiated energy 
which we estimate to be 3.408% ± 0.102% of the total 
ADM mass. The code is thus suitable for detailed stud- 
ies of various types of multiple black-hole simulation with 
regard to the generation of accurate waveform templates. 

In preparation for the comparison of black-hole col- 
lisions using different types of initial data, we test the 
code's performance and estimate in detail the error mar- 
gins associated with the different evolutions. Specifi- 
cally, we separately demonstrate convergence of the code 
for simulations starting from Brill-Lindquist, superposed 
Kerr-Schild and Misner data. We also study in depth 
the dependence of the resulting waveforms on the ex- 
traction radii. While the resulting uncertainties are rela- 
tively small for Brill-Lindquist and Misner data, we find 
the use of finite extraction radii to be the dominant error 
source for simulations of Kerr-Schild data. In the case 
of Brill-Lindquist data we further demonstrate the codes 
reliability by comparing head-on collisions with results 
available in the literature. 

We use the code to provide a detailed comparison of 
black-hole-binary head-on collisions using all three data 
types. In addition to the total mass of the system, ei- 
ther initial configuration has one free parameter which is 
specified by fixing the binding energy E-^/M of the sys- 
tem. We have compared the resulting waveforms for four 
initial configurations. 

The resulting waveforms obtained from Brill-Lindquist 
and Misner data show excellent agreement and predicts 
an energy radiated in the infall and merger of about 
0.052% - 0.056% of the total ADM mass M with the 
exact value increasing with the initial black-hole separa- 
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tion. In parallel the amount of energy due to spurious 
radiation inherent to the conformally flat initial data de- 
creases from 0.0031 % M to 0.0009 % M as the initial 
separation of the holes is increased from a binding en- 
ergy = -0.029 M to -0.0169 M. The uncertainties 
in these results for both data types are of the order of 
a few percent. While this good agreement might be ex- 
pected, given the similar nature of two data types, it is 
reassuring to confirm this expectation with high accuracy 
using current numerical techniques. 

In the case of the superposed Kerr-Schild data, we ob- 
serve a substantially larger amount of energy in the spuri- 
ous gravitational radiation due to the initial data, about 
a factor 15 larger than for the other data types. For 
sufficiently large separations, most of this spurious wave 
content radiates away before the merger of the holes and 
can thus be distinguished from the actual signal of the 
merger and ring-down. For the smaller separations stud- 
ied, however, this distinction becomes more difficult and 
leads to a non-ncgligiblo uncertainty in the amount of 
energy radiated in the infall and merger of the holes. 

In comparison with the conformally flat data types, the 
merger waveforms extracted at finite radius show larger 
amplitudes in the case of the Kerr-Schild simulations by 
5—10 %. The agreement of the resulting radiated energies 
in the infall and merger becomes much better, though, 
after extrapolating results to infinite extraction radius. 
Still, for large black- hole separations there remains a dis- 
crepancy of a few % in the merger energy between Kerr- 
Schild data on the one side and Brill-Lindquist and Mis- 
ner data on the other, even when taking into account 
remaining uncertainties in the simulations. We can there- 
fore not rule out systematic errors aflPecting the accuracy 
of the Kerr-Schild simulations. 

Such systematic errors can arise from the fact that the 
initial data are not inherently time symmetric and might 
imply small initial boosts of the individual holes. By 
evolving puncture data with nonvanishing Bowen-York 
momentum we have shown that small boosts can account 
for the discrepancies of the observed magnitude. A sec- 
ond systematic error arises from the constraint violations 
of the Kerr-Schild data. In particular, we observe the en- 
ergy contained in the spurious initial radiation, likely to 
be a signature of the constraint violations, to be signifi- 
cantly larger than the differences we observe. 

Finally, we mention future directions of research en- 
couraged by this study. First, it will be valuable to un- 
derstand the origin for the relatively large uncertainties 



in the wave amplitudes obtained from Kerr-Schild data 
at finite extraction radius or, conversely, why wave ex- 
traction appears to work remarkably well for the Misner 
and puncture simulations at radii significantly smaller 
than the wave zone. Second, it will be important to pro- 
duce evolutions of the constraint-solved version of the 
Kerr-Schild data and compare the results with those of 
the present study. An important question in this regard 
also concerns the amount of spurious initial radiation. 
A key advantage of Kerr-Schild-type initial data over 
conformally flat data is the fact that they contain the 
Kerr solution as a limit. They are therefore particularly 
promising candidates for producing initial data sets con- 
taining black holes with very large spins with minimal 
artificial radiation (see, for example, [91] for a discussion 
of artificial radiation in black-hole spacetimes with large 
spins). It is thus important to study how the solving 
of the constraints reduces the large amounts of artificial 
radiation observed here in the nonspinning case. 
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APPENDIX A: THE BSSN EVOLUTION 
EQUATIONS 

The BSSN equations are implemented in the Lean 

CODE using cither the set of variables defined in Eq. (1) 
or using the variable x i^i place of <j) as defined in Eq. (2). 
The (j) version of the BSSN equations used in Lean is 
given by 



dtiij = rdmlij + - -^lijdmr - '^ocA,^. (A1) 

dt4> = rdmct> + lidmf^ - aK), (A2) 

dtA,j = rdmA,, + 2i„(,9,)/3" - + e"^'^ (ai?^ - D.Djof^ + a{K i,, - 2i,'"i„,) , (A3) 
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dtK = rOruK - D^^D^a + a + -K^J , (A4) 

-\af^dmK + 2i'™ (6aa™</) - a„a) - ('^ + ^) (f ' - 7™"rL) S,/3^ (A5) 

I 



Here Di is the covariant derivative operator and Rij the 
Ricci tensor associated with the physical three-metric 7^ 
and the superscript denotes the trace- free part. We 
also note that the last term in the evolution equation 
(A5) vanishes in the continuum limit by virtue of the 
definition of in Eq. (1). With the addition of this term 
we follow Yo et. al. [92] who introduced this modification 



to improve the stability of the BSSN formulation in cases 
of relaxed symmetry assumptions of the spacetime under 
study. We set the free parameter a in this term to 2/3. 

The X version of the evolution system is obtained by 
substituting x = e~^'^. The evolution equations for the 
variables x, ^'■^ and F* are then given by 



dtX = 



dtAij 



2 T 



--af^dmK - A'-' 



3a 



dmX 



2d„ 



(A6) 
(A7) 

(A8) 



r 



while Eqs. (Al) and (A4) remain valid without modifi- 
cation. 



APPENDIX B: THE ADM VARIABLES OF A 
SINGLE BOOSTED KERR-SCHILD BLACK 
HOLE 

Purpose of this section is to calculate the ADM func- 
tions 7,j, Kij, a and /3* as functions of the laboratory 
coordinates for a nonspinning boosted black hole with 
mass parameter m in Kerr-Schild coordinates moving 
with speed in the laboratory rest frame^. The rest 
frame coordinates of the black hole are related to the 
laboratory coordinates by a Lorentz transformation 



where the transformation matrix is given by 



(Bl) 



(B2) 



^ Sec, for example, [93] for a discussion of the superposition of 
Kerr-Schild holes with nonvanishing spin. 



In the black- hole rest frame, the spacetime metric is given 
by (see e. g. [70]) 



9a0 = Voc0 + '^Hls.l-p, 



where 



r 



f = 



1, 



f 

X Xa-) 



(B3) 

(B4) 

(B5) 
(B6) 



and indices of and u° are raised and lowered with the 
flat space metric S^-^. 

The spacetime metric in the laboratory frame is ob- 
tained from that in the black-hole frame by a Lorentz 
transformation 

g^^ = A"^A^.g^0 = r?^. + 2He^£^, (B7) 

where we have used the fact that H and behave like a 
scalar and vector, respectively, and the Minkowski metric 
is invariant under Lorentz transformations. 

From the spacetime metric we directly obtain the 3- 
metric, its inverse as well as lapse and shift 

'Jmn — ^mn ~t~ '^-^^m^n^ (-^^) 

7"" = (5""-2i?5™'=(5"'4^i/[l-h2if(£o)^], (B9) 
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a = [l + 2//(4)^]-^/^ 

Z?'" = 2i7^o(5™'=4/[l + 2iJ(4)2] 



(BIO) 
(Bll) 
(B12) 



The extrinsic curvature is obtained from the derivatives 
of the three-metric according to 

= -^{dt7mn-2D^mPn)), (B13) 

where Cp denotes the Lie derivative along the shift vec- 
tor and -Dto the three-dimensional covariant derivative 
operator. 

The' dc^rivativos of the thrcc-mctric are most conve- 
niently expressed in terms of H and im 



dtlmn 



■Hi^dtin), (B14) 



2[V(™5„)ff- 



Finally, the derivatives of H and are given by 



= -A%^ 



(B15) 



(B16) 



b x A \a \b ^aXb^ (B17) 



In summary, the function for calculating the ADM vari- 
ables of a boosted black hole in Kerr-Schild coordinates 
at a particular point requires as input the coordinates of 
the point in the laboratory frame as well as the velocity 
V of the hole. First, the coordinates arc transformed into 
the rest frame of the black hole according to Eq. (Bl). 
Next H and follow from Eqs. (B4), (B5) and give the 
spacetime metric components (B7). The three- metric, 
lapse and shift follow from Eqs. (B8), (BIO) and (B12). 
Together with the extrinsic curvature (B13), these are 
returned to the calling function as the ADM variables of 
the boosted Kerr-Schild metric. 



APPENDIX C: THE ELECTROMAGNETIC 
DECOMPOSITION OF THE WEYL TENSOR 
AND WAVE EXTRACTION 

We calculate the Newman-Penrose scalar from the 
Weyl tensor via 



(CI) 



where n and m form part of a null-tetrad i, n, m, m such 
that all their inner products vanish except 



—i-n=l = m-m 



(C2) 



Specifically, we construct i, n and m from the orthonor- 
mal triad vectors u, v and w according to 

v2 



(C3) 



where fi^ is the timelike orthonormal vector. The triad 
u, V, w is constructed via Gram-Schmidt orthonormal- 
ization starting with 



u = [x ,y ,z\, 

v' = [xz , yz , -x'^ - y^] 



(C4) 



where e'™" represents the three-dimensional Levi-Civita 
tensor. 

In the decomposition of the Weyl tensor we follow the 

presentation of Friedrich [94] . The electric and magnetic 
part of the gravitational field are given by 



(C5) 
(C6) 



where -L^a = 5^a+n^na is the projector onto the space- 
like hypersurface and the * denotes the Hodge dual. By 
virtue of the Gauss-Codazzi equations (see, e. g. [95]), 
one can express the electromagnetic parts in terms of 
"3-1-1" variables according to 



Ty ^ kmn 7-) T/- 

J->ij — lik^ ^m-^nj • 



(C7) 



The Weyl tensor is then given in terms of the electric and 
magnetic part by Eq. (3.10) of Ref. [94]. Inserting this 
relation together with Eqs. (C3), (C7) into the definition 
(CI) enables us to express ^"4 exclusively in terms of 
"3-1-1" quantities 



^4=2 {Emn{v"'v^ - W"'W^ 



In practice, ^4 is calculated using this relation on the 

entire Cartesian grid and then interpolated onto coordi- 
nate spheres of different extraction radii. We then apply 
a mode decomposition using spherical harmonics of 
spin- weight —2 (cf. Eq. (42) in Ref. [96]) according to 



(C8) 



^A{t,e,4>)Y^:;^{e,<i>)d^i. (C9) 



In this context we note that 'I'4 is always extracted onto 
the entire coordinate sphere 9 — 0...7r, (j) = 0... 27r, 
even when underlying symmetry of the physical problem 
is used to reduce the computational domain to a bitant 
or octant. In those cases we use the fact that the real 
part of ^4 behaves like a scalar while the imaginary part 
of ^'4 behaves like a pscudo scalar, i. e. reverses its sign 
across symmetry boundaries. 
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TABLE IV: Performaiice summary for representative simula- 
tions presented in this work. 



Simulation 


dt/dx 


#CPU 


Memory 
[Gb] 


Speed 
[M/CPUh] 


Rl {h = 1/40) 


1/2 


12 


20.1 


0.620 


Rl {h = 1/44) 


1/2 


16 


26.4 


0.347 


Rl {h = 1/48) 


1/2 


24 


35.2 


0.225 


BL4 {h = 1/40) 


1/2 


8 


9.9 


0.854 


BL4 (/i = 1/44) 


1/2 


12 


12.5 


0.582 


BL4 {h = 1 /48) 


1/2 


16 


18.2 


0.350 


KS4 (ft = 1/20) 


1/4 


8 


8.1 


0.571 


KS4 (ft = 1/24) 


1/4 


8 


14.0 


0.315 


KS4 (ft = 1/28) 


1/4 


16 


22.5 


0.171 


M4 (ft = 1/320) 


1/8 


8 


9.1 


0.269 


M4 (ft= 1/360) 


1/8 


12 


12.3 


0.162 


M4 (ft = 1/400) 


1/8 


12 


17.5 


0.111 



APPENDIX D: PERFORMANCE OF THE CODE 

The majority of simulations presented in this work 
have been performed using a 24 node Linux cluster. Each 
node contains four AMD 2200 GHz processors and pro- 



vides 8 Gb of memory. Parallelization is implemented 
using the MPICH version 1.2.6 2004/08/08 libraries. The 
code has been compiled with version 4.0.2 2005091 of 
the gee, g + + and gf ortran compilers. Compared with 
alternative architectures (cf. Ref. [97]), we have noticed 
that this architecture requires about 25 % more memory 
resources for identical simulations. 

In Table IV we summarize the performance of the code 
for simulations Rl, BL4, KS4 and M4 of Table II using 
this architecture. The columns show the Courant factor 
dt/dx which scales linearly with the code's speed, the 
number of processors used in the simulation, the required 
memory as well as the speed. The latter is measured in 
physical time in units of the ADM mass M of the system 
per real time and processor. 

Regarding the memory usage we observe minor varia- 
tions, typically below 5 %, in the course of the simulation. 
This is due to the merger of refinement components as 
the black holes approach each other. The merger of re- 
finement components also leads to an increase in speed 
because the costly regridding operation is no longer re- 
quired and the total number of grid points of the merged 
refinement component is smaller than the sum of the two 
individual ones prior to merging. All reported speeds are 
averages over the entire simulation. 
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